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SUMMARY 


A compressible stability analysis computer code is 
developed. The code uses a matrix finite-difference method 
for local eigenvalue solution when a good guess for the 
eigenvalue is available and is significantly more compu- 
tationally efficient than the commonly used initial-value 
approach. The local eigenvalue search procedure also re- 
sults in eigenfunctions and, at little extra work, group 
velocities. A globally convergent eigenvalue procedure is 
also developed which may be used when no guess for the 
eigenvalue is available. The global problem is formulated 
in such a way that no unstable spurious modes appear so 
that the method is suitable for use in a black-box stability 
code. Sample stability calculations are presented for the 
boundary layer profiles of an LFC swept wing. 



SECTION 1 


INTRODUCTION 


The stability properties of compressible laminar boundary 
layers are particularly relevant to the phenomenon of laminar- 
turbulent flow transition. Recently, interest in this problem 
has increased because of applications to Laminar Flow Control 
(LFC) technology. In such applications there is a need for 
fast computer codes to perform efficient design calculations. 
The computer code SALLY [1] was developed for this purpose. 

It can perform optimized stability calculations for determin- 
ing suction requirements for maintaining laminar flow over 
swept wings. However, SALLY uses incompressible stability 
theory so it solves the eigenvalue problem for the fourth- 
order Orr-Sommerfeld differential equation. 

The linear stability analysis of three-dimensional com- 
pressible boundary layers involves solution of an eigenvalue 
problem for an eighth order system of differential equations. 
In the case of two-dimensional boundary layers or in the 
absence of dissipation in three-dimensional flow, the eighth- 
order system reduces to sixth order. 

The basic equations for the linear stability analysis of 
parallel-flow compressible boundary layers are derived using 
the samll disturbance theory. Infinitesimal disturbances of 
sinusoidal form are imposed on the steady boundary layer flow 
and substituded in the compressilbe Navier-Stokes equations. 
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Assuming that the mean flow is locally parallel, a set of 
five ordinary differential equations is obtained. Of these, 
there are three second order momentiom equations , one second 
order energy equation and one first order continuity equation. 
Following Lin [2-4] (his work involved only sixth order 
system) , this system is reduced to a set of eight first 
order ordinary differential equations making the system 
amenable to initial-value numerical integration procedures. 

All previous numerical investigations of compressible 
flow stability [5-10] make use of the initial value approach 
for the solution of this system of eight first-order equa- 
tions (or the reduced system of six first-order equations) . 

In these studies, the integration is started at a boundary 
(usually the free stream boundary since the asymptotic solu- 
tion of the stability equations provides initial values for 
starting the integration) and marched to the other boundary 
(solid wall) typically using a Runge-Kutta integration pro- 
cedure. Four linearly independent solutions are sought 
by means of this integration. The difficulty encountered 
in this scheme is that a straight-forward integration fails 
to produce four independent solutions because of parasitic 
error growth. This difficulty is usually overcome by an 
orthonormalization technique [see, e.g. , 11]. Upon obtaining 
four accurate linearly independent solutions, linear combina- 
tions are formed to satisfy all but one boundary condition 
at the wall (the boundary towards which the equations are 
marched) . The remaining boundary condition is satisfied only 
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when the eigenvalue condition is satisfied. It is normally 
imposed iteratively using a Newton-Raphson procedure which 
results in the desired eigenvalue. 

These initial-value methods to solve the compressible 
stability equations are often computationally slow and thus 
inefficient if used in a black-box stability code. Further- 
more, the initial-value methods require a good guess for the 
eigenvalue which is usually not available when one encounters 
a new problem. In such cases one can obtain eigenvalues using 
a local search procedure only by trial and error, which is 
neither elegant nor efficient. 

In the present work a computer code is developed for the 
compressible stability analysis of three-dimensional boundary 
layers. The code includes two eigenvalue search procedures— 
global (which is to be used when no guess is available) and 
local (which is used when a good guess is available) . The 
stability equations are solved in their original (3 second - 
order momentum equations, one second- order energy and one 
first-order continuity equation) using a matrix finite- 
difference' method. The reduction of the normal momentum 
equation to a first order equation for pressure fluctuations 
as done by Lees and Lin [2] is not done here since that re- 
sults in unstable spurious modes when the problem is solved 
using the global method (see below) . 
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The finite-difference representation of the stability 
equations results in a block- tridiagonal system of equations. 

A generalized matrix eigenvalue problem is then set up and 
solved using the complex LR algorithm [12] for global eigen- 
value search. The local search is performed by inverting the 
block-tridiagonal system using a block LU factorization to- 
gether with inverse Rayleigh iteration for the eigenvalues 
[12] in which the eigenvalue, eigenfunction and its adjoint 
are obtained simultaneously. One can obtain the group ve- 
locities which are needed in an eigenvalue optimization pro- 
cedure [13] at little extra cost to the local eigenvalue 
search. The procedure used in the present study for local 
eigenvalue solution is significantly faster than the initial 
value approach employed by previous investigators. 

Some sample calculations are performed for the stability 
analysis of compressible three-dimensional boundary layer 
profiles obtained for a laminar flow control wing using the 
boundary layer code developed by Kaups & Cebeci (see ref. [1]). 
particular wing chosen is a 35^^ swept super-critical wing 
with an 8 foot chord. The free stream Mach number is 0.89. 
Mack [7] and Lekoudis [9] reported stability calculations 
for the same wing using their stability codes. We have 
chosen three stations on the wing which represent forward 
and rearward crossflow instability regions and a midchord 
streamwise instability region. Some of the relevant boundary 
data is given in Table 1. More details on the wing and press- 
ure and suction distributions are given in reference [7] and 
[9]. 5 


SECTION 2 - COMPRESSIBLE STABILITY EQUATIONS 


Consider the stability of a three-dimensional locally- 
parallel compressible boundary layer flow. Let us use a 
cartesian coordinate system x,y,z in which the y-axis 
is normal to the solid boundary and x,z are parallel to 
it. (In the particular case of a wing, the x-axis will be 
assumed to be in the direction of the normal chord and the 
z-axis will be along the span.) If u,v,w are the x,y,z 
components of the instantaneous velocity, respectively, and 
p and T are instantaneous pressure and temperature, 
then, assuming that the base flow is locally parallel. 


u(x,y , z ,t) 
v(x,y ,z ,t) 
w(x,y ,x,t) 
p(x,y ,z,t) 
T (x,y ,z ,t) 


U^(y) + u(y)e 


i (ax+3z-wt) 


= v(y) e 


i (ax+3z-cot) 


(y) + w(y) e 


i (ax+3z-cot) 


= p(y) e 


i (ax+3z-o)t) 


T^Cy) + T(y)e 


i (ax+3z-wt) 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


Here represent the steady unperturbed boundary 

layer and quantities with tildas denote complex disturbance 
amplitudes . 

Substituting Eqs. (1-5) into the compressible Navier- 
Stokes equations, it can be shown that the linear parallel 
disturbances satisfy the following system of ordinary differ 
ential equations: 
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0 


( 6 ) 


(A + B D + C)J = 
where (j) is a five-element vector defined by 

{au + 3w, V, p, T, aw - 3u} (7) 

Here A,B,C are 5x5 matrices, and r where y is the 

normal boundary layer coordinate. The non- zero elements of 
the matrices A,B,C are given in Appendix I. 

The boundary conditions for Eq. (6) are 

y = 0; = <)>5 = 0 (8) 

<\>2f ^ 4 ^ 0 

Equations (6-8) constitute an eigenvalue problem for 
the frequency to and wavenumbers a, 3 . For a given Reynolds 
number R, this eigenvalue problem provides a complex dispersion 
relation of the form 

(0 = to (a , 3) (9) 

relating the parameters a, 3 and to which, in general, are all 
complex. In our analysis we choose to use temporal stability 
theory in which a, 3 are real and to is complex. It is thus 
assumed that the wavelike disturbances have x and z com- 
ponents of wave number a and 3 , respectively, and have a 
frequency to^ = Re (to) . It is further assumed that the dis- 
turbances grow or decay only in time. They grow if 
to^ = Im(co) >0 and decay if to^ < 0. 
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The selection of temporal theory results in a linear 
eigenvalue problem for o) so that global eigenvalue techniques 
can be applied. It is possible to extend the global method 
to the spatial problem for the compressible stability equa- 
tions but the work involved and the computer resources re- 
quired discourage any such attempt. It appears that the most 
efficient way to automatically provide a guess for the local 
solution of the spatial stability problem is to first solve 
the temporal global problem and then obtain the spatial guess 
using the group velocity transformations [14]. 

Equation (6) represents three second order momentum 
equations, one second-order energy equation and one first 
order continuity equation. It is possible to eliminate 
pressure entirely between these equations and to reduce 
the system to four second-order equations for velocity and 
temperature fluctuations. This reduction may seem desirable 
for the numerical integrations since the order of the matrix 
will be reduced. However, the elimination of pressure makes 
the equations singular and their solution becomes difficult 
to obtain in a numerically stable way. Similarly, it is 
possible to reduce the system of eight first-order equations 
to a single eighth-order differential equation [15], but 
this too is computationally ill conditioned. Thus, the follow- 
ing analysis is carried out in terms of five equations: three 

second-order equations for streamwise and spanwise velocity and 
for energy, one first-order equation for continuity and one 
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second order (but effectively first-order) equation for 
the normal velocity component. 



SECTION 3 


FINITE-DIFFERENCE REPRESENTATION 


The governing system of equations (6) for compressible 
flow is represented using a second order finite-difference 
formula on a staggered mesh (see Figure 1) . 

First the boundary layer coordinate Y , 0 < y < y is 
mapped into the computational domain 0 £ ^ £ 1 by the alge- 
braic mapping 


n 


g y 

L+y 


(10) 


where g = 1 + — 


Here y^ is the location of the edge of the boundary layer 

and L is a scaling parameter chosen to optimize the accuracy 

of the calculations. After some experience it has been found 

that a good choice for L is L 2y^ where y^ is the value 

of y at which the streamwise component of the mean velocity, 

U^, achieves 1/2 its freestream value. In the present study 

we choose y = 100. However, a much smaller value of y (=15) 
■^e “^e 

can be chosen for the local method since the free stream 
boundary conditions are imposed using the asymptotic solution. 

The computational domain n is divided into N equal inter- 
vals and the second order equations are represented as 


f^ A. 

1 3 


An 




1 


2 An 


+ c. 

3 
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] = 0 


+ d. 


[fsB. 


^ 2 ^2 

4- r* 

^2 ^ 2 

An j 

■ a 

J 

2 J 


(j=l,...,N - 1) (11) 


where (p . is the value of $ at n 
D 

(p, ■ (k = 1,...,5). Also, 

■1^ f D 

= 1, = 0 except 

d^ = 0 , d 2 = 1 for the p component 
of $, 


4 = = (g~n) 

1 2 2 
g L 


= j/N and has components 


f = 2 (g-n) 

2 „ 2_2 

g L 


f = (9-n) 

3 gL 


The first order continuity equation is represented as 

. + C . . 1 - n 

An :+TT = 0 


f„ B., 

3 D+2 


ii 


(12) 


( j - 0,. . . ,N - 1) 

Equations (11,12) along with the 8 boundary conditions 
(8) represent 5N+4 equations for 5N+4 unknowns. Since the 
velocity and temperature disturbances are assumed to be 
identically zero at the solid boundary (n = 0) the system 
reduces to 5N equations for 5N unknowns when these boundary 
conditions are applied. This is a block-tridiagonal system 
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of equations with 5x5 blocks which is solved using fully 
pivoted LU factorization. 

In the global method, the free-stream conditions (n = 1) 
are that cj)j^=(j)2==«J>^=<f>^=0 at n = !• This results in a linear 
eigenvalue problem for n. In the local method, asymptotic 
behavior ofJasy^“is found from (6) (see appendix II) . 
This asymptotic behavior is used to obtain a free-stream 
boundary condition of the form 

(ED + F) $ = 0 (13) 

that is applied at n = 1 on components cj)j^ (k=l , 2 , 4 ,5) . 
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SECTION 4 


GLOBAL METHOD 


When no guess is available for the eigenvalue of interest, 
it is best to use a method that is globally convergent and nearly 
guaranteed to converge to the eigenvalue. Such a method may 
be based on algorithms for calculations of the eigenvalues 
of a general complex matrix [12]. 

When the compressible stability equations (11-12) are for- 
mulated as a matrix problem, they take the form 

A $ = to B $ (14) 

where to is the eigenvalue and $ is the discrete representation 
of the eigenfunction. The eigenvalue is determined by the 
determinant condition 

Det| A - to B| =0 (15) 

For the present problem, B is invertible so (15) may be 
solved as 

Det I b”^ a - to 1 1 = 0 (16) 

which is the standard matrix eigenvalue problem, solvable by 
LR methods [ 12 ] . 

One has to be very careful in formulating the eigenvalue 
problem (14) using global methods to avoid the generation of 
growing (unstable) modes that are not physically relevant. 

These spurious unstable modes do not correspond to solutions 
of the differential equation— as the spatial resolution used 
to discretize the eigenfunction increases true modes of the 
differential equation converge while spurious modes do not. 
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A clumsy way to distinguish spurious modes from true 
modes is to change the spatial resolution and retain only 
those modes that do not change appreciably. This is neither 
efficient nor elegant. 

A better way is to eliminate the spurious unstable modes 
entirely. Spurious stable modes are still possible, but since 
these stable modes are normally very stable, they are not of 
much interest and can be easily disregarded without testing 
their true nature. We shall now describe a technique for 
eliminating the spurious unstable modes. 

The idea is simply that the spurious unstable modes 
would, if we used the same numerical method used for the 
stability problem on an initial-value problem instead, cause 
the unconditional instability of the numerical solution of 
the initial-value problem. On the other hand, if we were 
careful enough to use a numerical method for the stability 
problem that was also numerically stable for the initial- 
value problem, then no spurious unstable modes would exist. 

Thus, one way to avoid spurious modes is to set up the 
problem as one would for an initial value problem and to use 
a finite-difference scheme for the eigenvalue analysis that 
is consistent with the scheme for the initial-value problem. 
This is the method used in the present study; no spurious 
unstable modes appear in our calculations. In some initial 
work on this problem, we followed Lees and Lin [2] by re- 
ducing the second-order normal momentum equation to a first- 
order equation for pressure using the continuity equation. 
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This gives three second-order equations (two momentum and 
one energy equation) and two first-order equations (one for 
pressure and another for normal velocity) , When the matrix 
eigenvalue problem is set up using this system of equations 
by means of finite-differencing, several unstable spurious 
modes appear. 

When the eigenvalue problem (16) is solved using the 

2 

LR algorithm, the storage requirements are of 0(K ) while 

3 

the computational work involved is of 0 (K ) where K = 5N for 
the eighth order system of equations and K = 4N for the sixth- 
order system. It can be seen that the global method is 
relatively expensive computationally so it should be only 
used when no guess of the eigenvalue is available. 

Some computer timings for the global method on a CYBER 
175 computer are given in Table 2. All timings were obtained 
using the internal clock and are averaged over the three test 
cases listed in Table 1. Since the eigenvalue obtained by 
the sixth-order system is not much different from that of 
the eighth-order system, the use of the sixth-order system 
for the global method is recommended. 
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SECTION 5 - LOCAL EIGENVALUE SEARCH 


If a guess for the eigenvalue is available then one 
can improve its value by a local method. One way to perform 
the local analysis is to use a simple iterative method to 
find the eigenvalues of the matrix equation (15) that approxi- 
mates the compressible stability equations. An effective and 
efficient procedure for doing this is to use the inverse 
Rayleigh iteration procedure [12]. The generalization of 
this procedure to the compressible stability problem results 
in the following algorithm 


(A - B) J = B $ 

(A - B)"^ J (k+1) ^ ^ j (k) 

^ktl - ^j(k-M)^- J(ktl)j 


(17) 

(18) 
(19) 


The iteration cycle is started with a guessed eigenvalue 
and by assuming any smooth functional distribution for the 
eigenfunction J(0) and its adjoint J(0). In an integration 
of stability characteristics across a wing J(0) and if(0) may 
be chosen as the eigenfunction and adjoint from the previous 
station. 

At the end of each iteration cycle k the eigenfunction 
and its adjoint are normalized so that 
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( 20 ) 


The algorithm (17-20) has a rapid (cubic) rate of con- 
vergence. The error satisfies ^ ~ 0 ( (Wj^-to) ^) . 

VJe solve the block-tridiagonal form of Eq. (17) by 
using a fully pivoted LU method, in which case the same LU 
factorization applies to the solution of the adjoint problem 
(18) . 

In practice it is not necessary (rather, not recommended) 
to update the eigenvalue approximation (ji3j^ after each 
iteration to the eigenfunction and its adjoint using (17-18) . 
We have found it to be most efficient to iterate (17-18) 
approximately 4-10 times while keeping fixed (and, 

therefore, using the same LU factorization of A-to^^B) . Once 
the eigenfunction is refined less than 5 iterations are 
required for a fixed cOj^. Generally, only two outer itera- 
tions [LU factorizations and applications of (19) to update 
to required to converge to an eigenvalue. Some 

data on the speed of the local eigenvalue solution is given 
in Table 3. The time reported also includes the calcula- 
tion of group velocity. Since the eigenfunction and its 
adjoint are available as a result of local eigenvalue search, 
it costs little extra work to compute group velocity (see be- 
low) . 
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To show the sensitivity of the local method to the 
initial guess some sample calculations are made and re- 
ported in Table 4. The converged eigenfunction distribu- 
tions vs y for case 1 using 80 points are plotted in 
Figure 2 . 
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SECTION 6 - CALCULATION OF GROUP VELOCITY 


The group velocity is of importance in relating the 
results of spatial and temporal stability theory and in 
several optimization problems [13]. In a layered flow 
with three-dimensional disturbances having wavevector 
(a, 3) and frequency w(a,6) , the group velocity v is 


V 

g 


' 33 ^ 


( 21 ) 


One way to compute the group velocity is simply to 
compute the frequency 03 for several nearby values of a, 3 
and then use finite-difference approximations to v^. How 
ever, this procedure is not very efficient. 

A much better way to compute the group velocity is 
to first write the compressible stability equations for 
three dimensional disturbances in the form 


L(a, 3,w (a, 3) ) J = 0 


( 22 ) 


Taking the derivative of (22) with respect to cl gives 


3L 

9ot 




(23) 


Taking the inner product of (23) with the adjoint ^ of the 
eigenfunction $ gives 


3co _ (il/'3a 


(24) 
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since (J,L. /9ot) = d^/dct) = 0. There is a similar 

expression for Note that the inner product in (24) is 

the usual vector inner product; 

N 

(f,g) = ■ I (25) 

1=1 

It is not necessary to use either the adjoint eigenfunction 
of the differential equations (6) or the inner product 
/f (y) g (y) dy ; J is the adjoint of J with respect to the 
simple discrete inner product (25) and suffices to annihilate 
the term 90/ 9 a in (23) . 

The computation of the group velocity using (24) re- 
quires only the calculation of 9L/9a $(and 9L/93 0) since 

9L/9o) $ are available from the local eigenvalue itera- 

tions . 

In Table 5 , group velocities calculated by the present 
method for the three test cases are compared with those cal- 
culated by central differencing: 


9(jo I 

1 OJ . , T -03 . _ 

_ 1+1 1-1 

"9a 

Ij '"j+l-'^j-l 

9(jo 


9B 

Ij 


which required four extra eigenvalue calculations so the 
cost of determining group velocity was four times the cost 
of the eigenvalue search. Using the present method, group 
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I 


velocity can be obtained at less than 10% of the cost of the 
eigenvalue search. 
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SECTION 7 


RICHARDSON EXTRAPOLATION 


The finite difference method presented above is only 
second order accurate. However, the accuracy of the eigen- 
value (and the group velocity) obtained by this procedure 
can be enchanced by Richardson extrapolation to the limit 
An = 1/N ->0 [15] . 


If 


S^?^ = w(h.) (i = 0,...m) 


(27) 


where h^ is the grid size and w is the eigenvalue computed 

2 

by a method with error of 0(h ) , then 


q(j-l) _ s^3-l) 

+ ^-1 i 

1 1+1 h. 9 

‘hTTT) - 1 
1+3 


(28) 


(j=l,....,m, i=0, #in-j) 

The resulting approximation has an error of order 

0(h^'^^) when appropriate sequences of grid sizes h^ are used. 

We present here eigenvalue calculations for different grids 

(h. = ^ ) ^nd their extrapolated values. The sequence of 
i 

h^ chosen is that proposed by Bulirsch and Stoer [16]. 

Tables 6,7,8 give eigenvalue results for boundary layer 
data cases 1,2 and 3, respectively. 

It can be seen that an eigenvalue converged to 5 sig- 
nificant digits is obtained as a result of the extrapolation 
procedure. For most applications, an eigenvalue that has 
converged to 3 significant digits should suffice. This is 
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achieved by extrapolating only between three eigenvalues 
at N = 20, 40 and 60. For three point extrapolation, (28) 
reduces to the following formula 

^extrapolated = | [h^-h^) /h^ ]W q + 

[{h^-h2)/h2]Uj^ + [(h^-h^)/h^)002}/ 

I [(h^-h2)/hp] + [(h2-h2)/h^] 

+ [ (hg-h^)/h2] I (29) 

It is our experience that a reasonable answer can be 
obtained by this extrapolation formula using only 20, 30 
and 40 points. Eigenvalues obtained using this set of 
points are compared with the converged value using 20, 40, 
60, 80, 120 and 160 points in Table 9. The eigenvalue 
accurate to 3 significant digits is obtained in about 2 
seconds of CYBER 175 time. 


SECTION 8 - CONCLUSION 

Efficient local and global eigenvalue methods for the 
stability analysis of plane-parallel three-dimensional com- 
pressible boundary layers have been developed. The computer 
code that implements these matrix boundary-value methods is 
significantly more efficient than previous codes based on 
initial-value methods and is suitable for a black-box sta- 
bility analysis package. 
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APPENDIX I 


The Stability Equation 



The stability equations for three-dimensional compressible 
boundary layers are written as in (6) : 

(A + B D + C) $ = 0 

where 

cj)j^ = a u + 3 w 

^2 = V 

<t'3 = P 

*4 = T 

(p^ = a w - 3 u 

and 



Here u, v, w represent the complex amplitudes of the per- 
turbation velocity in x, y, z directions, respectively and p, 

T are the perturbation amplitudes of pressure and temperature. 

The disturbances are assumed to be of the form 
J(y) + gz - (ot) 

where a, 3 are x and z wave numbers, respectively and o) is the 
(complex) frequency. 

The non-zero elements of 5x5 matrices A, B and C are given 
by 


A 

A 

A 

A 


11 

22 

44 

55 


1 

1 

1 

1 


26 



I dy 

^11 y dT o 
®12 ^ 


-i- ioLV' + 3W") 

y dT o o 

o O 


®21 " 

B = i ^ 
22 dT^ 

R 


®32 = 1 


2(y- 1)M^ a (au^ + 3W^)/(a + 3 ) 

o dvi_ 


B = ^ ^ T 

44 y^ dT^ o 


= 2(y- 1)M^ a (aW^ - 3U^)/(ct + B ) 


B - i ^ 

54 “ dT^ 

B = -i- ^ 

55 dT^ 

Q = .. ^ 

11 ^y T 


r 

^12 " ■“ *■ y T~ 


(aW^ - BU^) 


2 2 

(aU + 3W - 03) + X (a + B ) 1 
o o 

^ -• • '" 2 . 2 

(aU^ + BW^) - ^ ^ (a +6)1 

(a^ +6^) 


o 


, 1 

,2 

d 

(aUo + 

3W 
, o 

' ^ 

dT^ 

o 

. X-2 

dy 

o 

T 



dT^ 

o 

o 


_ r ^ 

L . , rp 

R 

\ 

(aU 

+ 3W 

f 
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^ dy ^ -* 

I f dr; ew^) 

O O 


^31 ^ 


32 


C^2 = i Y M - 0)) 

C^= - T~ (ctU + 3W^ - w) 

34 T o o 

o 


42 


- [5_g__ - 2 i (y-I)M^ a(aU^ + 3W^) ] 


y T o 
o o 


'4 3 


i-5-2: (y-i)M^ (au^ + 3W^ - (O) 


C 44 — [ 


^ ^ (aU + 3W - 0)) + (a^ + 3^) 


y T 

o 


o 


o , dy 

- (Y-i) o M rr df“ ) 

o 


2 ^2 
1 

o o 


2 1 


,2 

T Cl y ^ ^ ^ m 

— — (T ) - U„ dT„ o ] 

y ^nn 2 O ^00 

O dT ^ 
o 


'52 


— ^ (aw^ - 3U^) 

y T o o 


■54 y 


o o 

1 O m 

2 o 


o dT 


. " 1 *^^0 

(aW - su^) + (“«o 

o o 


3U^ ) 


^55 ~ ^ 


i R 

o o 


(aU + 3W - 0 }) + a + 3 ) ] 


Here subscript o refers to the unperturbed boundary layer 
and primed quantities represent the derivative of the quantity 
with respect to the boundary layer coordinate y. 


28 



In the above, y and y are the dynamic viscosity and 
the ratio of the specific heats respectively for the gas which 
is assumed to be perfect. The Prandtl number a is assumed to 
be constant. Moreover, X is defined as 

X = I + 2) 


where \X2 is the ratio of the second coefficient of viscosity 
to the first. 

In the present study all velocities have been scaled 
by U^, the x component of velocity at the edge of the boundary 

•k 

layer, and all lengths are scaled by 6 , the displacement thick- 
ness of the velocity profile in x-direction, U^(y). The result- 
ant Reynolds number and Mach number are then given by 
* 



e 


U 

M = ^ 

/y A T 

e 

where and T^ are the kinematic viscosity and mean temperature 
in the free stream, and A is the universal gas constant. 

The results reported in the present paper were obtained 
with a = 0.72 and y^ = 1.2. Dynamic viscosity, y, was cal- 
culated using Sutherland's viscosity law. At least for the 
test cases considered in the present work, the results are 
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not sensitive to the value of y 2 * 

The special choice of = au + 3w and <f)^ = aw - 3u 
allows the eighth-order system of equations to reduce to 
the sixth order for eigenvalue calculation if is assumed 

to be zero. This assumption implies the absence of dissipa- 
tion. 
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APPENDIX II 


BOUNDARY CONDITIONS 



In the free stream (y ^ governing equations (6) 


reduce to 


D (j)^ + + C,, (j)^ = 0 


11 '^1 '^13 ^3 


D^(f)^ + D 4 >t + D (f)^ + (J)^ = 0 


23 


22 ^2 


2 '"31 ^1 ^33 ^3 

D"c|>y, + 4^4 = 0 


D<t2 + C3j^ (l)^ + C33 4.3 + C34 <f>4 = 0 
2 


D 4)5 + C53 (j>3 = 0 


(II-l) 

(II-2) 

(II-3) 

(II-4) 

(II-5) 


J 

where D = #-. The non-zero coefficients B. . and C. . are the same 
dy 13 13 

as in appendix I except that they are no more dependent on y. By 


eliminating 4>^, the above equations can be written as 
D^i|j_ + a Dip^ + eijj- + = 0 




where 


and 


D 1 P 4 + q4^4 = 0 


4^1 = (pj_, 4^2 = 4 ^ 2 ' 4^3 = ^4' ^4 = 4^5 


(II- 6 ) 

(II-7) 

(II- 8 ) 

(II-9) 


^ “ ®12 ^13^^33 


b = (B 


21 


^23^31 

^33 


B 


23, 


)/(l - 

^33 


32 



B 


C /{I 

^ " "^ 43-^^33 


® ^11 ^13 ^ 31^^^33 


^ ^13 ^ 34^^33 


g 

h 


022/(1-323/033 ) 

^43 ^31 


s = O 


44 


33 
- 0 


43 ^34^^^33 


q = 0 


55 


Equations (II- 6 ) - (II-9) admit a general solution of 
the form 


8 (j) Xj(y-y ) 

4 - = E p z e ^ ® ;i=l, 

j=i :j ^ 


(II-IO) 


where X . are the eight characteristic roots of the coefficient 
matrix of Equations (II- 6 ) - (II-9) and are the components 

of the eight characteristic vectors 
The asymptotic condition 
0 as y 00 

is imposed by requiring that the arbitrary constants p^ 
vanish when 


Re(A^) > 0. 
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There are four which satisfy this condition resulting 
in four boundary conditions. 

The roots (j = can be obtained from the 

characteristic equation (derived from Eqs. (II-6) - (II-8)) 

6 4 2 

X + L 2 X + = 0 

where 

= g + e + k - ab ” cd 

= ke - hf + eg + kg - abk + bdf - ced + ahc 
= keg - ghf 

2 

Assuming that X = X , the above sixth order equation 
reduces to the cubic equation 

+ L 2 X + = 0 


The roots of this cubic equation can be found as in 
[17] and consequently six of the X^'s are known. The com- 
ponents of the corresponding characteristic vectors are 


(j) 


(j) 


1 

d X? + de - ah 
3 

ka - fd + aX^ 

3 


= - (X? + e + f z^^^)/a X. ; j = 1,6 
2 2 3 J 


The seventh and eighth characteristic root is simply 
given as 

X^ = (-q) ^ 

Xg = -(-q)*' 
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and the corresponding components of the characteristic 
vectors are 


z<’>= 1 
4 

1 


All other components z are zero. 

i 

Now, let 


M = 


K 1 


where 


E p . X . 2 : 


(j) 


(y-y^) 




j=l 


j j 1 


;i=l,4 


so 


8 (j) >>j(y-y ) 

m. = Z n,. ' e ^ ® ;i=l;8 

j=l D 1 


where 


N = 


z<3> 

1 


X .Z ^ ^ ^ 

L ) 


;i-l,4 


The arbitrary constants p^ are obtained by inverting 
8x8 matrix N, 

P = M 

The asymptotic boundary conditions require that P vanish 
whenever Re(Xj)>0. This results in four equations of the 
form 

(ED + F) ip = 0 

where both E and F are 4x4 matrices. 
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Table 1.- Characteristics of the sample test cases on 

a 35^ swept wing 




(chord. 

c = 8 ft) 

• 



Case 

No. 

Streamwise 

Location, 

x/c 

R 

M 

a 

3 

t 

1 

0.001868 

145 

0.386 

0.272117 

-0.29181 

51.96° 

2 

0.4639 

3615 

1.058 

0.1153 

0.0 

26.93° 

3 

0.8921 

1754 

0.736 

0.2432 

-0.2654 

34.80° 

t Angle 

formed by the 

local potential flow vector with 

the x-axis. 




Table 2,- Timings for the global eigenvalue 
method (time given in seconds on a 
CYBER 175 computer) . 



8th order 

6th order 

N 

system 

system 

15 

3.15 

2.05 

20 

7.12 

5.17 

25 

13.47 

8.65 


Table 3 

Timings for the local 
method 

eigenvalue 


8th order 

6th order 

N 

system 

system 

20 

0.61 

0.40 

40 

1.14 

0.77 

80 

2.22 

1.52 
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Table 4. -Effect of the initial guess on the 
convergence of the local eigenvalue method 
for case 3 using N = 80, (8th-order system). 


Initial guess 

Time required 
for convergence 

(0.011595,0.002 3522) * 

1.82 

(0.01,0.002) 

2.38 

(0.02,0.004) 

2.95 

(0.02,0.008) 

3.10 

(0.005,0.001) 

3.11 

(0.03,0.006) 

3.50 

(0.04,0.004) 

not converged 


* equals the converged value 



Table 5. 
cases 1-3 

-Computed group velocities for 
using N = 80, (8th-order system) 

. 


bi 

a 

"^3 


Case 

No. 

Central 

difference 

formula 

Present 

method 

Central 

difference 

formula 

Present 

method 

1 

(0.66039, 

-0.00625) 

(0.66036, 

-0.00636) 

(0.73769 

0.01386) 

(0.73771, 

0.01398) 

2 

(0.39893, 

0.02105) 

(0.39876, 

0.02063) 

(0.21884, 

-0.000987) 

(0.21867, 

-0.00155) 

3 

(0.53507, 

-0.05212) 

(0.53512, 

-0.05219) 

(0.44301, 

-0.04472) 

(0.44295, 

-0.04464) 


40 





Table 6. -Richardson extrapolation of the most unstable eigenvalue for case 1 


i 

N m= 

0 

1 

Re{w) X 10 
2 

3 

4 

5 

1 

20 

-.2686117 

-.2629801 

-.2630158 

-.2630098 

-.2630087 

-.2630100 

2 

40 

-.2643880 

-.2630118 

-.2630102 

-.2630087 

-.2630099 


3 

60 

-2.636234 

-.2630106 

-.2630089 

0.2630099 



4 

80 

-.2633553 

-.2630093 

-.2630097 




5 

120 

-.2631631 

-.2630096 





6 

160 

-.2630959 







i 

N m= 

0 

1 

Im((jo) X 100 
2 

3 

4 

5 

1 

20 

.6073290 

.6191717 

.6186855 

.6185509 

.6186694 

.6186523 

2 

40 

.6162110 

.6187395 

.6185593 

.6186661 

.6186525 


3 

60 

.6176158 

.6186043 

.6186542 

.6186534 



4 

80 

.6180483 

.6186418 

.6186535 




5 

120 

.6183780 

.6186506 







6 


160 


6184972 



Table 

7. - Richardson 

extrapolation 

of the most 

unstable eigenvalue 

for case 2 


i 

N 

m= 0 

1 

Re((jo) X 10 
2 

3 

4 

5 

1 

20 

.3940779 

.3877680 

. 3889774 

. 3888798 

.3888478 

.3888521 

2 

40 

.3893454 

.3888430 

. 3888859 

.3888486 

.3888520 


3 

60 

.3890663 

.3888751 

.3888528 

.3888518 



4 

80 

.3889827 

.3888584 

.3888519 




5 

120 

. 3889136 

. 3888535 





6 

160 

.3888873 


Im{(jo) X 100 




i 

N 

o 

1 

2 

3 

4 

5 

1 

20 

.0725901 

.1230929 

.1228928 

.1227987 

.1228500 

.1228501 

2 

40 

.1104672 

.1229150 

.1228046 

.1228486 

.1228501 


3 

60 

.1173826 

.1228322 

.1228437 

.1228500 



4 

80 

.1197668 

.1228408 

.1228491 




5 

120 

.1214746 

.1228470 






6 


160 


1220750 



Table 8, -Richardson extrapolation of the most unstable eigenvalue for case 3 






Re ( (Ji3 ) X 10 




i 

N m= 

0 

1 

2 

3 

4 

5 

1 

20 

.1171532 

.1158833 

.1158746 

.1158764 

.1158752 

.1158759 

2 

40 

.1162008 

.1158756 

.1158763 

.1158752 

.1158759 


3 

60 

.1160201 

.1158761 

.1158754 

.1158758 



4 

80 

.1159571 

.1158756 

.1158758 




5 

120 

.1159118 

.1158757 





6 

160 

.1158960 


Im(co) X 100 




i 

N m= 

0 

1 

2 

3 

4 

5 

1 

20 

.2200906 

.2362895 

.2362199 

.2362203 

.2362156 

.2362186 

2 

40 

.2322398 

.2362276 

.2362203 

.2362158 

.2362185 


3 

60 

.2344553 

.2362221 

.2362163 

.2362184 



4 

80 

.2352283 

.2362177 

.2362181 




5 

120 

.2357780 

.2362180 





6 

160 

.2359705 








Table 9 .-Comparison of eigenvalues obtained 
using different orders of extrapolation 


Case 

No. 


Three point 
extrapolation 
using N = 20, 
30 and 40 


Six point 
extrapolation 
(see tables 
6 - 8 ) 


1 

(-0.026303, 

(-0.026301, 


0.0061848) 

0.0061865) 

2 

( 0.038889, 

( 0.038885, 


0.0012324) 

0.0012285) 

3 

( 0.011587, 

( 0.011587, 


0.0023623) 

0.0023621) 



( FREE STREAM ) 


77=1 


1/2 

I 


0 


(m 



-I nr.v.w.Y) 

-j + l/2 
— j + I 


N-l 

N-l/2 


77 =0 
(WALL) 


Figure 1.- A SCHEMATIC OF THE STAGGARD FINITE-DIFFERENCE GRID 
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